Annual model - Regions
Tanya Rogers
12/15/2021
library(dplyr)
#library(psych) #for pairs.panels, but could use other packages, e.g. GGalley
library(lavaan)
library(semPlot)
library(DiagrammeR)
library(tidyr)
library(ggplot2)Import data
combined=read.csv("data/annual_averages/annual_data_compiled_regions.csv",stringsAsFactors = F)
cnames=read.csv("analysis/column_names_region.csv", stringsAsFactors = F)
dsub=filter(combined, Year>=1975) %>% arrange(Region,Year)
focaldata=dsub[,cnames$Datacolumn]
fvars=cnames$Shortname
colnames(focaldata)=fvars
regions=unique(focaldata$region)
regionorder=c("West","North","South")
years=1975:2021
focaldata$tzoop=focaldata$hcope+focaldata$clad+focaldata$mysid+focaldata$pcope
fvars=c(fvars,"tzoop")
cnames=rbind(cnames,data.frame(Longname=NA,Shortname="tzoop",
Diagramname="Total Zooplankton\nBiomass",
Datacolumn=NA,Log="yes"))
#focal variables
varnames=c("temp", "flow","nitrate","ammonia","dophos","chla","hcope","clad","amphi","pcope","mysid","potam","corbic","sside","estfish_bsot","estfish_bsmt","tzoop")
source("analysis/myLavaanPlot.r")Data prep
Log transform, scale
#log transform
logvars=fvars[cnames$Log=="yes"]
logtrans=function(x) {
x2=x[which(!is.na(x))]
if(any(x2==0)) {log(x+min(x2[which(x2>0)],na.rm=T))}
else {log(x)}
}
focaldatalog = focaldata %>%
mutate_at(logvars,logtrans)
#scale data
fdr0=focaldatalog
tvars=fvars[-(1:2)]
fdr=fdr0 %>% group_by(region) %>%
#lag
mutate_at(tvars,list("1"=lag)) %>%
#scale
mutate_at(-(1:2),scale) %>%
ungroup() %>%
as.data.frame()
#detrended data
fdr_dtr=fdr0 %>% group_by(region) %>%
#detrend
mutate_at(tvars,function(x) {
x<<-x
if(!all(is.na(x))) {
if((length(which(x==0))/length(x))<0.5) {
x2<<-x
x2[x2==0]=NA
res<<-residuals(lm(x2~years))
out=x
out[which(!is.na(x2))]=res
return(out)
} else {return(x)}
} else {return(x)}
}) %>%
#lag
mutate_at(tvars,list("1"=lag)) %>%
#scale
mutate_at(-(1:2),scale) %>%
ungroup() %>%
as.data.frame() Time series plots
Cross-correlation matrices
(only sig correlations shown… no correction for multiple comparisons)
Fish indices are not correlated in S and N!
SEM model
With and without detrending.
#west has no ssides, corbic
# modwest='zoop=~hcope+mysid
# fish=~estfish_bsmt+estfish_bsot
# zoop~chla+potam+flow
# chla~potam+flow
# fish~zoop+flow
# '
modwest='chla~potam+flow
tzoop~chla+potam+flow
estfish_bsmt~tzoop+flow
estfish_bsot~tzoop+flow
'
modfitwest=sem(modwest, data=filter(fdr,region=="West"))
modfitwest_dtr=sem(modwest, data=filter(fdr_dtr,region=="West"))
summary(modfitwest, standardized=T, rsq=T)## lavaan 0.6-4 ended normally after 27 iterations
##
## Optimization method NLMINB
## Number of free parameters 14
##
## Used Total
## Number of observations 40 47
##
## Estimator ML
## Model Fit Test Statistic 4.813
## Degrees of freedom 4
## P-value (Chi-square) 0.307
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## potam -0.325 0.135 -2.407 0.016 -0.325 -0.380
## flow 0.098 0.131 0.748 0.454 0.098 0.118
## tzoop ~
## chla 0.831 0.096 8.676 0.000 0.831 0.792
## potam -0.154 0.088 -1.756 0.079 -0.154 -0.171
## flow -0.069 0.080 -0.867 0.386 -0.069 -0.080
## estfish_bsmt ~
## tzoop 0.626 0.115 5.457 0.000 0.626 0.568
## flow 0.386 0.100 3.874 0.000 0.386 0.403
## estfish_bsot ~
## tzoop 0.598 0.119 5.030 0.000 0.598 0.542
## flow 0.389 0.103 3.771 0.000 0.389 0.407
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish_bsmt ~~
## .estfish_bsot 0.054 0.066 0.811 0.418 0.054 0.129
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .chla 0.584 0.131 4.472 0.000 0.584 0.802
## .tzoop 0.215 0.048 4.472 0.000 0.215 0.267
## .estfish_bsmt 0.402 0.090 4.472 0.000 0.402 0.412
## .estfish_bsot 0.430 0.096 4.472 0.000 0.430 0.442
##
## R-Square:
## Estimate
## chla 0.198
## tzoop 0.733
## estfish_bsmt 0.588
## estfish_bsot 0.558summary(modfitwest_dtr, standardized=T, rsq=T)## lavaan 0.6-4 ended normally after 24 iterations
##
## Optimization method NLMINB
## Number of free parameters 14
##
## Used Total
## Number of observations 40 47
##
## Estimator ML
## Model Fit Test Statistic 6.849
## Degrees of freedom 4
## P-value (Chi-square) 0.144
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## potam 0.032 0.162 0.196 0.844 0.032 0.034
## flow 0.215 0.160 1.341 0.180 0.215 0.229
## tzoop ~
## chla 0.837 0.105 8.009 0.000 0.837 0.795
## potam 0.026 0.107 0.247 0.805 0.026 0.026
## flow -0.020 0.108 -0.183 0.854 -0.020 -0.020
## estfish_bsmt ~
## tzoop 0.417 0.115 3.618 0.000 0.417 0.428
## flow 0.446 0.114 3.915 0.000 0.446 0.463
## estfish_bsot ~
## tzoop 0.379 0.121 3.139 0.002 0.379 0.383
## flow 0.455 0.119 3.824 0.000 0.455 0.467
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish_bsmt ~~
## .estfish_bsot 0.046 0.089 0.519 0.604 0.046 0.082
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .chla 0.879 0.197 4.472 0.000 0.879 0.953
## .tzoop 0.384 0.086 4.472 0.000 0.384 0.376
## .estfish_bsmt 0.534 0.119 4.472 0.000 0.534 0.548
## .estfish_bsot 0.584 0.131 4.472 0.000 0.584 0.585
##
## R-Square:
## Estimate
## chla 0.047
## tzoop 0.624
## estfish_bsmt 0.452
## estfish_bsot 0.415# par(mfrow=c(1,2))
# semPaths(modfitwest, "std", edge.label.cex = 1, residuals = F)
# semPaths(modfitwest, "par", edge.label.cex = 1, residuals = F)
labelswest <- createLabels(modfitwest, cnames)
# residuals(modfitwest)
# modificationindices(modfitwest)# modnorth='zoop=~hcope+mysid
# #fish=~estfish_bsmt+estfish_bsot
# zoop~chla+potam+flow
# chla~potam+flow
# estfish_bsmt~zoop+flow
# estfish_bsot~zoop+flow
# '
# modnorth='zoop=~clad
# zoop~chla+corbic+potam+flow
# chla~corbic+potam+flow
# estfish_bsmt~zoop+flow+sside+chla
# estfish_bsot~zoop+flow+sside+chla
# '
modnorth='chla~corbic+potam+flow
tzoop~chla+corbic+potam+flow
estfish_bsmt~tzoop+flow+chla
estfish_bsot~tzoop+flow
'
modfitnorth=sem(modnorth, data=filter(fdr,region=="North"))
modfitnorth_dtr=sem(modnorth, data=filter(fdr_dtr,region=="North"))
summary(modfitnorth, standardized=T, rsq=T)## lavaan 0.6-4 ended normally after 22 iterations
##
## Optimization method NLMINB
## Number of free parameters 17
##
## Used Total
## Number of observations 40 47
##
## Estimator ML
## Model Fit Test Statistic 1.857
## Degrees of freedom 5
## P-value (Chi-square) 0.868
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## corbic 0.358 0.132 2.718 0.007 0.358 0.422
## potam -0.121 0.143 -0.844 0.399 -0.121 -0.140
## flow 0.022 0.138 0.159 0.874 0.022 0.026
## tzoop ~
## chla 0.681 0.112 6.072 0.000 0.681 0.614
## corbic 0.414 0.102 4.069 0.000 0.414 0.440
## potam 0.075 0.103 0.735 0.463 0.075 0.078
## flow -0.351 0.098 -3.586 0.000 -0.351 -0.370
## estfish_bsmt ~
## tzoop 0.033 0.188 0.175 0.861 0.033 0.032
## flow 0.040 0.130 0.304 0.761 0.040 0.041
## chla 0.694 0.215 3.227 0.001 0.694 0.614
## estfish_bsot ~
## tzoop 0.468 0.135 3.476 0.001 0.468 0.470
## flow 0.251 0.128 1.968 0.049 0.251 0.266
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish_bsmt ~~
## .estfish_bsot 0.025 0.098 0.256 0.798 0.025 0.040
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .chla 0.562 0.126 4.472 0.000 0.562 0.737
## .tzoop 0.283 0.063 4.472 0.000 0.283 0.301
## .estfish_bsmt 0.565 0.126 4.472 0.000 0.565 0.579
## .estfish_bsot 0.677 0.151 4.472 0.000 0.677 0.727
##
## R-Square:
## Estimate
## chla 0.263
## tzoop 0.699
## estfish_bsmt 0.421
## estfish_bsot 0.273summary(modfitnorth_dtr, standardized=T, rsq=T)## lavaan 0.6-4 ended normally after 23 iterations
##
## Optimization method NLMINB
## Number of free parameters 17
##
## Used Total
## Number of observations 40 47
##
## Estimator ML
## Model Fit Test Statistic 2.521
## Degrees of freedom 5
## P-value (Chi-square) 0.773
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## corbic 0.315 0.162 1.940 0.052 0.315 0.311
## potam 0.066 0.153 0.431 0.667 0.066 0.070
## flow 0.105 0.169 0.619 0.536 0.105 0.108
## tzoop ~
## chla 0.710 0.097 7.337 0.000 0.710 0.679
## corbic 0.317 0.104 3.049 0.002 0.317 0.299
## potam 0.086 0.094 0.920 0.358 0.086 0.087
## flow -0.460 0.104 -4.414 0.000 -0.460 -0.454
## estfish_bsmt ~
## tzoop 0.356 0.181 1.972 0.049 0.356 0.373
## flow 0.204 0.135 1.507 0.132 0.204 0.211
## chla 0.333 0.187 1.784 0.074 0.333 0.334
## estfish_bsot ~
## tzoop 0.228 0.149 1.528 0.126 0.228 0.238
## flow 0.225 0.151 1.492 0.136 0.225 0.232
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish_bsmt ~~
## .estfish_bsot 0.085 0.109 0.776 0.438 0.085 0.124
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .chla 0.846 0.189 4.472 0.000 0.846 0.871
## .tzoop 0.317 0.071 4.472 0.000 0.317 0.298
## .estfish_bsmt 0.526 0.118 4.472 0.000 0.526 0.542
## .estfish_bsot 0.892 0.199 4.472 0.000 0.892 0.916
##
## R-Square:
## Estimate
## chla 0.129
## tzoop 0.702
## estfish_bsmt 0.458
## estfish_bsot 0.084# par(mfrow=c(1,2))
# semPaths(modfitnorth, "std", edge.label.cex = 1, residuals = F)
# semPaths(modfitnorth, "par", edge.label.cex = 1, residuals = F)
labelsnorth <- createLabels(modfitnorth, cnames)
# residuals(modfitnorth)
# modificationindices(modfitnorth)#no potam
# modsouth='zoop=~hcope+mysid
# #fish=~estfish_bsmt+estfish_bsot
# zoop~chla+corbic+flow
# chla~corbic+flow
# estfish_bsmt~zoop+flow
# estfish_bsot~zoop+flow
# '
modsouth='chla~corbic+flow
tzoop~chla+corbic+flow
estfish_bsmt~tzoop+flow+corbic
estfish_bsot~tzoop+flow+corbic
'
modfitsouth=sem(modsouth, data=filter(fdr,region=="South"))
modfitsouth_dtr=sem(modsouth, data=filter(fdr_dtr,region=="South"))
summary(modfitsouth, standardized=T, rsq=T)## lavaan 0.6-4 ended normally after 20 iterations
##
## Optimization method NLMINB
## Number of free parameters 16
##
## Used Total
## Number of observations 40 47
##
## Estimator ML
## Model Fit Test Statistic 1.175
## Degrees of freedom 2
## P-value (Chi-square) 0.556
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## corbic 0.446 0.154 2.902 0.004 0.446 0.436
## flow -0.014 0.147 -0.093 0.926 -0.014 -0.014
## tzoop ~
## chla 0.677 0.095 7.107 0.000 0.677 0.705
## corbic 0.242 0.102 2.371 0.018 0.242 0.246
## flow -0.138 0.088 -1.559 0.119 -0.138 -0.147
## estfish_bsmt ~
## tzoop -0.008 0.164 -0.052 0.959 -0.008 -0.008
## flow -0.020 0.140 -0.142 0.887 -0.020 -0.021
## corbic 0.491 0.170 2.896 0.004 0.491 0.496
## estfish_bsot ~
## tzoop 0.184 0.150 1.231 0.218 0.184 0.191
## flow -0.175 0.128 -1.366 0.172 -0.175 -0.193
## corbic 0.423 0.155 2.727 0.006 0.423 0.446
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish_bsmt ~~
## .estfish_bsot 0.004 0.108 0.035 0.972 0.004 0.006
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .chla 0.846 0.189 4.472 0.000 0.846 0.813
## .tzoop 0.307 0.069 4.472 0.000 0.307 0.320
## .estfish_bsmt 0.745 0.167 4.472 0.000 0.745 0.765
## .estfish_bsot 0.623 0.139 4.472 0.000 0.623 0.697
##
## R-Square:
## Estimate
## chla 0.187
## tzoop 0.680
## estfish_bsmt 0.235
## estfish_bsot 0.303summary(modfitsouth_dtr, standardized=T, rsq=T)## lavaan 0.6-4 ended normally after 15 iterations
##
## Optimization method NLMINB
## Number of free parameters 16
##
## Used Total
## Number of observations 40 47
##
## Estimator ML
## Model Fit Test Statistic 0.825
## Degrees of freedom 2
## P-value (Chi-square) 0.662
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## corbic 0.160 0.164 0.973 0.330 0.160 0.157
## flow -0.032 0.158 -0.206 0.837 -0.032 -0.033
## tzoop ~
## chla 0.639 0.112 5.708 0.000 0.639 0.642
## corbic 0.200 0.117 1.707 0.088 0.200 0.198
## flow -0.188 0.112 -1.685 0.092 -0.188 -0.193
## estfish_bsmt ~
## tzoop -0.223 0.153 -1.460 0.144 -0.223 -0.229
## flow -0.084 0.149 -0.566 0.571 -0.084 -0.089
## corbic 0.343 0.158 2.168 0.030 0.343 0.347
## estfish_bsot ~
## tzoop 0.082 0.146 0.562 0.574 0.082 0.085
## flow -0.220 0.142 -1.549 0.121 -0.220 -0.236
## corbic 0.349 0.151 2.316 0.021 0.349 0.360
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish_bsmt ~~
## .estfish_bsot -0.056 0.130 -0.429 0.668 -0.056 -0.068
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .chla 1.010 0.226 4.472 0.000 1.010 0.977
## .tzoop 0.506 0.113 4.472 0.000 0.506 0.494
## .estfish_bsmt 0.857 0.192 4.472 0.000 0.857 0.879
## .estfish_bsot 0.781 0.175 4.472 0.000 0.781 0.829
##
## R-Square:
## Estimate
## chla 0.023
## tzoop 0.506
## estfish_bsmt 0.121
## estfish_bsot 0.171# par(mfrow=c(1,2))
# semPaths(modfitsouth, "std", edge.label.cex = 1, residuals = F)
# semPaths(modfitsouth, "par", edge.label.cex = 1, residuals = F)
labelssouth <- createLabels(modfitsouth, cnames)
# residuals(modfitsouth)
# modificationindices(modfitsouth)Nice plots
Without covariances
Original units
West
North
South
Detrended
West
North
South
With covariances
Original units
Detrended